library(funFEM)
library(reshape2)
library(ggplot2)
library(gridExtra)
library(FactoMineR)
library(factoextra)
library(corrplot)
library(mclust)
library(stringr)
library(cluster)
library(clusterSim)
library(ggmap)

1 Etude préliminaire

1.1 Préparation des données

Cet ensemble de données contient des données provenant du système de partage de vélos de Paris, appelé Vélib. Les données sont des profils de chargement des stations de vélos sur une semaine. Les données ont été collectées toutes les heures pendant la période du dimanche 1er septembre au dimanche 7 septembre 2014.

On charge ici les données Velib disponibles dans la librairie funFEM.

library(funFEM)
data(velib)

Ces données contiennent

  • data : les profils de chargement (nb de vélos disponibles / nb de points d’attache) des 1189 stations à 181 points de temps.

  • position : la longitude et la latitude des 1189 stations de vélos.

  • dates : les dates de téléchargement.

  • bonus : indique si la station est sur une colline (bonus = 1).

  • names : les noms des stations.

Afin d’avoir des semaines complètes, nous allons supprimer les 13 premières colonnes.

Velibdata <- velib$data[, -c(1:13)]
colnames(Velibdata) <- velib$dates[-c(1:13)]

On controle le nombre de jours et d’heures dans le jeu de données

day <- str_sub(colnames(Velibdata), 1, 3)
day <- factor(day, levels = c("Lun", "Mar", "Mer", "Jeu", "Ven", "Sam", "Dim"))
table(day)
day
Lun Mar Mer Jeu Ven Sam Dim 
 24  24  24  24  24  24  24 
hour <- str_sub(colnames(Velibdata), 5, 6)
table(hour)
hour
00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
 7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7 
timeTick <- 24 * (0:6)

On stocke les informations sur les stations dans le dataFrame suivant :

station <- data.frame(nom = velib$names, velib$position, colline = velib$bonus)

On ajoute des informations complémentaires sur les stations :

# Les gares
gare <- rep(0, nrow(station))
gare[which(str_detect(station$nom, "GARE") == T)] <- 1

# Lieux touristiques (autour de Tour eiffel, Panthéon, Arc de Triomphe, Champs
# Elysées, Louvre, Notre Dame, Invalides, Trocadéro, Sacré Coeur, Montmartre)
tourisme <- rep(0, nrow(station))
TI <- c(209, 504, 897, 102, 518, 357, 256, 301, 992, 627, 579, 241, 646, 13, 283,
    941, 976, 1189, 22, 633, 853, 973, 762, 195, 71, 131, 912, 302, 303, 1068, 774,
    384, 307, 326, 436, 789, 250, 755, 461, 439, 918, 1065, 99, 151, 313, 871)
tourisme[TI] <- 1

# Les mairies
mairie <- rep(0, nrow(station))
mairie[which(str_detect(station$nom, "MAIRIE") == T)] <- 1
mairie[which(str_detect(station$nom, "HOTEL") == T)] <- 1

# Bois / Parc / Square
parc <- rep(0, nrow(station))
parc[which(str_detect(station$nom, "BOIS DE") == T)] <- 1
parc[which(str_detect(station$nom, "PARC") == T)] <- 1
parc[which(str_detect(station$nom, "SQUARE") == T)] <- 1
parc[which(str_detect(station$nom, "CANAL") == T)] <- 1

# Porte
porte <- rep(0, nrow(station))
porte[which(str_detect(station$nom, "PORTE") == T)] <- 1

station <- data.frame(station, gare = gare, tourisme = tourisme, mairie = mairie,
    parc = parc, porte = porte)

1.2 Quelques fonctions auxiliaires pour exploiter les résultats

Fonction pour tracer les profils (charge de chaque station / loading) en fonction des jours et heures pour les stations choisies.

plotprofils <- function(Velibdata, station, numstation, plotsem = TRUE) {
    # Velibdata = les données de velib initiales station = ens. des
    # informations sur les stations numstation = vecteur contenant le numéro de
    # ligne des stations dont on souhaite tracer les profils

    library(reshape2)
    library(ggplot2)
    Dataaux <- data.frame(id.s = station$nom, Velibdata)
    Aux <- melt(Dataaux[numstation, ], id = c("id.s"))
    Aux <- data.frame(Aux, day = str_sub(Aux$variable, 1, 3), hour = as.numeric(str_sub(Aux$variable,
        5, 6)))
    Aux$day <- factor(Aux$day, levels = c("Lun", "Mer", "Mar", "Sam", "Jeu", "Dim",
        "Ven"))

    if (plotsem == TRUE) {
        g <- ggplot(Aux, aes(x = hour, y = value, col = id.s)) + geom_line() + facet_wrap(~day,
            ncol = 2) + xlab("Time") + ylab("Loading") + scale_x_continuous(breaks = seq(0,
            24, 4)) + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
            labs(col = "Classif.") + scale_x_continuous(breaks = c(0, 7, 9, 12, 14,
            20, 24))
    } else {
        v <- rep(0, nrow(Aux))
        for (j in 1:length(levels(Aux$day))) v[which(Aux$day == levels(Aux$day)[j])] <- Aux$hour[which(Aux$day ==
            levels(Aux$day)[j])] + (24 * (j - 1))
        Aux <- data.frame(Aux, v = v)
        rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(0, 7),
            ymax = rep(1, 7), color = rainbow(n = 7))
        g <- ggplot(data = Aux) + geom_line(data = Aux, aes(x = v, y = value, col = id.s)) +
            geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
                fill = color), show.legend = FALSE, alpha = 0.1) + xlab("Time") +
            ylab("Loading") + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
            labs(col = "Classif.") + scale_x_continuous(breaks = rep(seq(0, 21, 3),
            7) + rep(24 * (0:6), each = 4), labels = rep(seq(0, 21, 3), 7))
    }
    return(g)
}

Exemple d’utilisation :

plotprofils(Velibdata, station, numstation = c(813, 655, 468))

plotprofils(Velibdata, station, numstation = c(813, 655, 468), plotsem = FALSE)

Se donnant une classification, fonction pour tracer les profils moyens par classe :

plotprofilsmoy <- function(Velibdata, clustering, plotsem = TRUE) {
    library(reshape2)
    library(ggplot2)

    Dataaux <- matrix(0, nrow = max(clustering), ncol = ncol(Velibdata))
    for (k in 1:max(clustering)) {
        Dataaux[k, ] <- apply(Velibdata[which(clustering == k), ], 2, mean)
    }
    colnames(Dataaux) <- colnames(Velibdata)
    Dataaux <- data.frame(id.s = as.factor(1:max(clustering)), Dataaux)
    Aux <- melt(Dataaux, id = c("id.s"))
    Aux <- data.frame(Aux, day = str_sub(Aux$variable, 1, 3), hour = as.numeric(str_sub(Aux$variable,
        5, 6)))
    Aux$day <- factor(Aux$day, levels = c("Lun", "Mer", "Mar", "Sam", "Jeu", "Dim",
        "Ven"))

    if (plotsem == TRUE) {
        g <- ggplot(Aux, aes(x = hour, y = value, col = id.s)) + geom_line() + facet_wrap(~day,
            ncol = 2) + xlab("Time") + ylab("Loading") + ylim(0, 1) + theme(legend.position = "bottom") +
            theme(axis.title.x = element_blank()) + labs(col = "Classif.") + scale_x_continuous(breaks = c(0,
            7, 9, 12, 14, 20, 24))
    } else {
        v <- rep(0, nrow(Aux))
        for (j in 1:length(levels(Aux$day))) v[which(Aux$day == levels(Aux$day)[j])] <- Aux$hour[which(Aux$day ==
            levels(Aux$day)[j])] + (24 * (j - 1))
        Aux <- data.frame(Aux, v = v)
        rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(0, 7),
            ymax = rep(1, 7), color = rainbow(n = 7))
        g <- ggplot(data = Aux) + geom_line(data = Aux, aes(x = v, y = value, col = id.s)) +
            geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
                fill = color), show.legend = FALSE, alpha = 0.1) + xlab("Time") +
            ylab("Loading") + ylim(0, 1) + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
            labs(col = "Classif.") + scale_x_continuous(breaks = rep(seq(0, 21, 3),
            7) + rep(24 * (0:6), each = 4), labels = rep(seq(0, 21, 3), 7))
    }
    return(g)
}

Exemple d’utilisation avec les stations sur colline par rapport aux autres :

plotprofilsmoy(Velibdata, clustering = station$colline + 1, plotsem = FALSE)

Fonction auxiliaire pour comparer une classification avec les informations auxiliaires disponibles pour les stations

stationcaract <- function(station, clustering) {
    library(dplyr)
    # Dataaux<-data.frame(id.s=c(1:nrow(Data)),station) aux <- cbind(Dataaux,
    # clust)
    aux <- cbind(station, clustering)
    aux.long <- melt(data.frame(lapply(aux[, -c(2:3)], as.character)), stringsAsFactors = FALSE,
        id = c("nom", "clustering"), factorsAsStrings = T)
    # Effectifs
    aux.long.q <- aux.long %>%
        group_by(clustering, variable, value) %>%
        mutate(count = n_distinct(nom)) %>%
        distinct(clustering, variable, value, count)
    # avec fréquences
    aux.long.p <- aux.long.q %>%
        group_by(clustering, variable) %>%
        mutate(perc = count/sum(count)) %>%
        arrange(clustering)

    gaux <- ggplot(data = aux.long.p, aes(x = clustering, y = perc, fill = value)) +
        geom_bar(stat = "identity") + facet_wrap(~variable)

    return(gaux)
}

1.3 Quelques statistiques descriptives

On reprend ici rapidement quelques éléments d’analyse descriptive vus précédemment dans l’UF.

  • Tracé des corrélations des heures par jour, toutes les stations confondues :
library(corrplot)
timeTickAux <- c(timeTick, 168)
for (k in 1:7) corrplot(cor(Velibdata[, (timeTickAux[k] + 1):(timeTickAux[k + 1])]),
    method = "ellipse")

  • Comportement moyen par jour :
velibmeanday <- matrix(0, nrow = nrow(Velibdata), ncol = 7)
for (k in 1:7) {
    velibmeanday[, k] <- apply(Velibdata[, (timeTickAux[k] + 1):(timeTickAux[k +
        1])], 1, mean)
}
colnames(velibmeanday) <- c("Lun", "Mar", "Mer", "Jeu", "Ven", "Sam", "Dim")
ggplot(melt(velibmeanday), aes(x = Var2, y = value)) + geom_boxplot() + xlab("") +
    ylab("Loading")

  • Comportement moyen par heure et par jour :
velibHour <- data.frame(value = colMeans(Velibdata), day = day, hour = as.numeric(hour))
ggplot(velibHour, aes(x = hour, y = value, col = day)) + geom_line() + geom_point() +
    ylab("Loading") + ggtitle("Hourly loading, averaged over all stations")

  • Quelques visualisations des stations sur le plan 2D de Paris :
library(ggmap)
# Average loading for each station
load <- rowMeans(Velibdata)
qmplot(longitude, latitude, data = station, color = load) + scale_colour_gradient(high = "red",
    low = "yellow")

# scale_colour_gradient(high='blue',low='white')
# Loading moyen à l'heure t (le premier temps est 0h)
t <- 0
loadheure <- rowMeans(Velibdata[, seq(from = (t + 1), by = 24, length = 7)])

qmplot(longitude, latitude, data = station, color = loadheure) + scale_colour_gradient(high = "red",
    low = "yellow")

  • Visualisation des stations selon les informations auxiliaires sur les stations :
qmplot(longitude, latitude, data = station, color = as.factor(station$colline))

qmplot(longitude, latitude, data = station, color = as.factor(station$tourisme))

1.4 Analyse en composantes principales

Question : Faites une ACP des stations et analysez les résultats.

respca <- PCA(.......)
fviz_eig(....)
fviz_pca_var(......)
fviz_pca_ind(........)

2 Clustering sur les coordonnées de l’ACP

Dans cette partie, nous allons utiliser les coordonnées de l’ACP comme matrice de données pour classer les stations.

velibacp <- respca$ind$coord[, 1:7]

2.1 Clustering avec l’algorithme des Kmeans

Question : A l’aide d’une procédure des Kmeans, déterminez une classification des stations.

# A COMPLETER

Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue.

# A COMPLETER

2.2 Clustering avec de la classification hiérarchique

Question : A l’aide d’une procédure hiérarchique, déterminez une classification des stations.

# A COMPLETER

Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue par CAH. Comparez cette classification avec celle obtenue par les Kmeans.

# A COMPLETER

2.3 Clustering avec des mélanges gaussiens

Question : A l’aide des mélanges gaussiens, déterminez une classification des stations.

# A COMPLETER

Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue par modèles de mélange. Comparez cette classification avec celles obtenues précédemment.

# A COMPLETER

3 Clustering sur un résumé des données

Dans cette partie, on va s’intéresser à un clustering sur un jeu de données “résumé” au sens suivant : on fait la moyenne des loadings selon les tranches horaires suivantes : 0-7h, 8h-20h et 21h-23h pour chaque jour et chaque station.

3.1 Préparation des données

VelibMean <- matrix(0, nrow = nrow(Velibdata), ncol = 3 * 7)
for (jour in 1:7) {
    VelibMean[, 3 * (jour - 1) + 1] <- apply(Velibdata[, 24 * (jour - 1) + c(1:8)],
        1, mean)
    VelibMean[, 3 * (jour - 1) + 2] <- apply(Velibdata[, 24 * (jour - 1) + c(9:21)],
        1, mean)
    VelibMean[, 3 * (jour - 1) + 3] <- apply(Velibdata[, 24 * (jour - 1) + c(22:24)],
        1, mean)
}
colnames(VelibMean) <- paste(rep(unique(day), each = 3), "-", rep(c("0-7h", "8h-20h",
    "21-23h"), 7), sep = "")

Question : Etudiez les corrélations entre variables du jeu de données VelibMean.

# A COMPLETER

Question : Faites une ACP et commentez.

# A COMPLETER

3.2 Méthodes de classification

3.2.1 Kmeans

Question : Faites une classification des stations à partir du jeu de données VelibMean à l’aide des Kmeans. Etudiez la classification obtenue.

# A COMPLETER

3.2.2 Mélanges gaussiens

Question : Faites une classification des stations à partir du jeu de données VelibMean à l’aide des mélanges gaussiens. Etudiez la classification obtenue.

# A COMPLETER

4 Clustering pour des séries temporelles

Les données de Vélib sont des données fonctionnelles (on étudie un phénomène qui évolue au cours du temps). Il existe des méthodes de clustering dédiées aux données fonctionnelles, comme par exemple le package funFEM.

Pour utiliser cette méthode, on commence par considérer la décompostion dans la base de Fourier des courbes.

basis <- create.fourier.basis(c(0, 167), nbasis = 25)
fdobj <- smooth.basis(1:168, t(Velibdata), basis)$fd

La méthode funFEM est une procédure de clustering basée sur des mélanges fonctionnels. Comme pour les mélanges gaussiens, il y a différentes formes disponibles selon les contraintes imposés dans le modèle (le détail est hors programme !). Cette méthode est détaillée dans l’article de Bouveyron et al. (2015).

Question : Executez les commandes suivantes et exploitez les résultats. On considère une classification en 10 classes comme dans Bouveyron et al. (2015).

resfunFEM <- funFEM(fdobj, K = 10, model = "AkjBk", init = "kmeans", lambda = 0,
    disp = TRUE)

df <- data.frame(proba = apply(resfunFEM$P, 1, max), classif = as.factor(apply(resfunFEM$P,
    1, which.max)))
ggplot(df, aes(x = classif, y = proba)) + geom_boxplot()
  • Etude de la classification obtenue
classifFEM <- apply(resfunFEM$P, 1, which.max)
plotprofilsmoy(Velibdata, clustering = classifFEM, plotsem = FALSE)

qmplot(longitude, latitude, data = station, colour = as.factor(classifFEM))

stationcaract(station, classifFEM)
  • Comparaison avec les classifications précédentes
table(resICL$classification, classifFEM)
adjustedRandIndex(resICL$classification, classifFEM)
table(resICLbis$classification, classifFEM)
adjustedRandIndex(resICLbis$classification, classifFEM)